library(tidyverse)
library(sf)
library(sfnetworks)
library(tidygraph)
library(igraph)
library(lwgeom)
library(nngeo)
library(mapview)
library(osmdata)
library(randomcoloR)
library(viridis)
library(gifski)
library(fs)
Sys.setlocale(locale = "hebrew")
gush_dan <- st_read("D:/Downloads/israel-ISR_gpkg/tel_aviv-4309.gpkg",layer = "edges")
dat1 <- opq_osm_id (type = "relation", id = "1382494") %>%
opq_string () %>%
osmdata_sf ()
tlv <- dat1$osm_multipolygons
step1 <- st_intersection(gush_dan,tlv)
step2 <- step1 %>% st_transform(2039)
step3 <- step2 %>% select()
step4 <- step3 %>% st_cast("LINESTRING")
pts <- step4 %>%
st_line_sample(density = 1/300) %>%
st_cast("POINT") %>%
st_sf() %>%
filter(!st_is_empty(geometry)) %>%
`$`(geometry)
step5 <- step4 %>%
as_sfnetwork() %>%
st_network_blend(pts,0.1) %>%
activate(edges) %>%
st_as_sf()
net <- as_sfnetwork(step5)
in_giant <- net %>% components() %>% `$`(membership) %>% `==`(1)
net1 <- net %>%
filter(in_giant) %>%
activate(edges) %>%
filter(!edge_is_multiple()) %>%
filter(!edge_is_loop()) %>%
activate(nodes)
mat <- net1 %>%
st_network_cost()
K <- 20
cols <- randomColor(K)
colnames(mat) <- rownames(mat) <- as.numeric(1:nrow(mat))
vec <- mat[,1:K] %>% apply(1,which.min)
new_centers <- map_dbl(1:K,function(x){
if(length(mat[vec == x,vec == x]) == 1){
x
}else{
mat[vec == x,vec == x] %>%
apply(1,max) %>%
which.min() %>%
names() %>%
as.numeric()
}
})
old_centers <- vec
p <- ggplot() +
geom_sf(data = net1 %>% st_as_sf() %>% mutate(vec = factor(vec)),mapping = aes(color = vec,fill = vec),stroke = 0.3,size = 2,shape = 21) +
geom_sf(data = net1 %>% st_as_sf() %>% slice( new_centers),color = "black",size =3) +
ggtitle("1") +
scale_color_manual(values = cols,guide = "none")
ggsave("01.png")
Z <- 100
lst <- list(Z)
i = 2
for(x1 in 1:Z){
print(x1)
lst[[x1]] <- net1 %>%
activate(nodes) %>%
st_as_sf() %>%
mutate(vec = vec, z=x1)
# reassign to clusters
vec <- mat[,new_centers] %>% apply(1,which.min)
mins <- mat[,new_centers] %>% apply(1,min)
print(max(mins))
# refind new centers
new_centers <- map_dbl(1:K,function(x){
if(length(mat[vec == x,vec == x]) == 1){
new_centers[x]
}else{
mat[vec == x,vec == x] %>%
apply(1,max) %>%
which.min() %>%
names() %>%
as.numeric()
}
})
if(identical(old_centers,vec)){
break
}
old_centers <- vec
p <- ggplot() +
geom_sf(data = net1 %>% st_as_sf() %>% mutate(vec = factor(vec)),mapping = aes(color = vec,fill = mins),stroke = 0.3,size = 2,shape = 21) +
geom_sf(data = net1 %>% st_as_sf() %>% slice( new_centers),color = "black",size =3) +
ggtitle(i)+
scale_color_manual(values = cols,guide = "none")+
scale_fill_viridis(option = "magma",limits = c(0,5000))
ggsave(paste0(ifelse(i<10,paste0("0",i),i),".png"))
i <- i + 1
}
paths <- dir_info() %>%
filter(str_detect(path,"png")) %>%
pull(path)
gif_file <- file.path("result.gif")
gifski(paths[-1], gif_file,width = 2100,height = 2100,delay = 0.2)
unlink(paths)
